if !d.name eq 'PS' then begin
    device,xsize=18,ysize=24,yoffset=3
    !p.thick=1 & !x.thick=1 & !y.thick=1
end

@eu_read_bcprofl
;!p.multi=[0,3,2]
;!p.charsize=2.5
!p.multi=[0,3,2] & !p.charsize=2.5 & !x.margin=[3.75,0] & !y.margin=[3,0.5]

default,nn,5
ntimes = (size(u))[3]
good = (ntimes-nn)+indgen(nn)-1
print,'doing averages over',nn, ' xaver records ...'

ruu = (total((total(u2,1)/m),1)/l)-(total((total(u^2,1)/m),1)/l)
rvv = (total((total(v2,1)/m),1)/l)-(total((total(v^2,1)/m),1)/l)
rww = (total((total(w2,1)/m),1)/l)-(total((total(w^2,1)/m),1)/l)


uu=fltarr(l,m)
vv=fltarr(l,m)
ww=fltarr(l,m)
omega=fltarr(l,m)

quu=fltarr(l,m)
qvv=fltarr(l,m)
qww=fltarr(l,m)
qwv=fltarr(l,m)
qwu=fltarr(l,m)
qvu=fltarr(l,m)
rsinth=fltarr(l,m)
rurms2=fltarr(l,m)

urms2t = ruu+rvv+rww
urms2 = (total(urms2t[good]/float((size(good))[1])))

; reynolds stresses
for k=0, l-1 do begin
  for j=0,m-1 do begin
    rsinth[k,j]=r[k]*sin(theta[j])
    uu[k,j] = total(u[j,k,good])/float((size(good))[1])
    vv[k,j] = total(v[j,k,good])/float((size(good))[1])
    ww[k,j] = total(w[j,k,good])/float((size(good))[1])
    quu[k,j]= total(u2[j,k,good]- u[j,k,good]^2)/float((size(good))[1])
    qvv[k,j]= total(v2[j,k,good]- v[j,k,good]^2)/float((size(good))[1])
    qww[k,j]= total(w2[j,k,good]- w[j,k,good]^2)/float((size(good))[1])
    qwv[k,j]= total(rwv[j,k,good]- w[j,k,good]*v[j,k,good])/float((size(good))[1])
    qwu[k,j]=total(rwu[j,k,good]-  w[j,k,good]*u[j,k,good])/float((size(good))[1])
    qvu[k,j]=total(rvu[j,k,good]-  v[j,k,good]*u[j,k,good])/float((size(good))[1])
    rurms2[k,j]=quu[k,j]+qvv[k,j]+qww[k,j]
  endfor
endfor

omega[*,1:m-2]=uu[*,1:m-2]/rsinth[*,1:m-2]

nlev=64
thg=where(th GE -85 AND th LT 85)

lev = 2.*max(abs(omega[*,thg]))*indgen(nlev)/float(nlev-1)-max(abs(omega[*,thg]))
contour,smooth(omega[*,thg],2),x[*,thg],y[*,thg],/fi,nl=64,/iso,title='!7X!6',lev=lev
colorbar,range=[min(lev),max(lev)],pos=[0.15,0.7,0.16,0.85],/vert,ytickformat='(F6.1)',yticks=2,ytickv=[min(lev),0.,max(lev)],yaxis=0,char=2.5

lev = 2.*max(abs(vv[*,thg]))*indgen(nlev)/float(nlev-1)-max(abs(vv[*,thg]))
contour,smooth(vv[*,thg],2),x[*,thg],y[*,thg],/fi,nl=64,/iso,title='!8v!6',lev=lev
colorbar,range=[min(lev),max(lev)],pos=[0.485,0.7,0.495,0.85],/vert,ytickformat='(F6.2)',yticks=2,ytickv=[min(lev),0.,max(lev)],yaxis=0,char=2.5

lev = 2.*max(abs(ww[*,thg]))*indgen(nlev)/float(nlev-1)-max(abs(ww[*,thg]))
contour,smooth(ww[*,thg],2),x[*,thg],y[*,thg],/fi,nl=64,/iso,title='!8w!6',lev=lev
colorbar,range=[min(lev),max(lev)],pos=[0.82,0.7,0.83,0.85],/vert,ytickformat='(F6.3)',yticks=2,ytickv=[min(lev),0.,max(lev)],yaxis=0,char=2.5

;lev = grange(-0.02,0.02,64)
lev = 2.*max(abs(qwv/rurms2))*indgen(nlev)/float(nlev-1)-max(abs(qwv/rurms2))
contour,smooth(qwv[*,thg]/rurms2[*,thg],2),x[*,thg],y[*,thg],/fi,nl=64,/iso,title='!8Q!Dr!7h!N!6',levels=lev
colorbar,range=[min(lev),max(lev)],pos=[0.15,0.2,0.16,0.35],/vert,ytickformat='(F6.2)',yticks=2,ytickv=[min(lev),0.,max(lev)],yaxis=0,char=2.5

lev = 2.*max(abs(qwu/rurms2))*indgen(nlev)/float(nlev-1)-max(abs(qwu/rurms2))
contour,smooth(qwu[*,thg]/rurms2[*,thg],2),x[*,thg],y[*,thg],/fi,nl=64,/iso,title='!8Q!Dr!7u!N!6',levels=lev
colorbar,range=[min(lev),max(lev)],pos=[0.485,0.2,0.495,0.35],/vert,ytickformat='(F6.2)',yticks=2,ytickv=[min(lev),0.,max(lev)],yaxis=0,char=2.5

lev = 2.*max(abs(qvu/rurms2))*indgen(nlev)/float(nlev-1)-max(abs(qvu/rurms2))
contour,smooth(qvu[*,thg]/rurms2[*,thg],2),x[*,thg],y[*,thg],/fi,nl=64,/iso,title='!8Q!D!7hu!N!6',levels=lev
colorbar,range=[min(lev),max(lev)],pos=[0.82,0.2,0.83,0.35],/vert,ytickformat='(F6.2)',yticks=2,ytickv=[min(lev),0.,max(lev)],yaxis=0,char=2.5

end
